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Abstract 

We consider a class of multiscale parabolic problems with diffusion coefficients oscillating in space at a possibly 
small scale e. Numerical homogenization methods are popular for such problems, because they capture efficiently 
the asymptotic behaviour as e —> 0, without using a dramatically fine spatial discretization at the scale of the 
fast oscillations. However, known such homogenization schemes are in general not accurate for both the highly 
oscillatory regime e —> 0 and the non oscillatory regime e ~ 1. In this paper, we introduce an Asymptotic 
Preserving method based on an exact micro-macro decomposition of the solution which remains consistent for 
both regimes. 

Resume 

Schemas numeriques Asymptotic Preserving pour les problemes paraboliques multi-echelles. On 

considere une classe de problemes paraboliques multi-echelles dont les coefficients de diffusion oscillent rapidement 
en espace a une echelle e possiblement petite. Les methodes numeriques d’homogeneisation sont populaires pour 
ces problemes, car elles capturent efficacement le comportement asymptotique lorsque e — > 0, sans utiliser une 
discretisation spatiale aussi fine que l’echelle des oscillations rapides, comme le necessiteraient les methodes non- 
raides standards. Cependant, les schemas d’homogeneisation existants ne sont en general pas precis dans les deux 
regimes oscillant e —> 0 et non-oscillant e ~ 1. Dans ce travail, nous introduisons une methode Asymptotic 
Preserving basee sur une decomposition micro-macro exacte qui reste consistante pour les deux regimes. 
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Version frangaise abregee 


L’objectif est de construire des schemas numeriques pour des problemes paraboliques (1) ou les coeffi¬ 
cients de diffusion sont hautement oscillants en espace. Ce type de rnodele intervient dans des problemes 
de propagation en milieux exhibants une structure periodique. Dans de nombreux cas, la taille de la 
periode est petite par rapport a la taille caracteristique du milieu, et en notant par e leur rapport, une 
analyse asymptotique est necessaire pour etudier le comportement de la solution quand e —> 0. Cette 
analyse a ete conduite dans [4,5,3,10] a l’aide de la convergence double-echelle. 

D’un point de vue numerique, une approche directe s’avere tres coiiteuse puisque les parametres 
numeriques des maillages doivent resoudre la plus petite echelle e. Les methodes numeriques d’ho- 
mogeneisation comrne la “heterogeneous multiscale method” (HMM) [7] (voir [2] dans le context de 
problemes paraboliques lineaires, et Particle de revue [1]) permettent de calculer efficacement la solution 
homogeneisee u° mais aussi la solution oscillante u £ a e fixe en se basant sur l’approximation du rnodele 
asymptotique, qui suppose que e est tres petit, ainsi que des techniques de correcteurs [5,10]. Cependant, 
si e n’est pas petit, ce type d’approche tombe en defaut. 

Dans ce travail, nous proposons une methode qui permet d’approclier la solution u £ pour n’importe 
quelle valeur de e £ (0,1], a parametres numeriques fixes independament de e. Ce type d’approche est 
dit ’’Asymptotic Preserving” [11] : un tel schema est consistant avec le probleme initial pour tout e fixe 
et degenere quand e —> 0 en un schema numerique consistant avec le rnodele asymptotique. 

Notre approche est basee sur une reformulation du probleme initial en un probleme augmente (voir 
(5) et les notations (6)) satisfait par U £ (t,x,y ), dans lequel les echelles lentes x et rapides y = x/e sont 
considerees independantes. La solution u e (t,x) du probleme initial peut alors etre retrouvee grace a la 
relation U £ (t,x,y = x/e) = u £ (t,x). Dans le probleme (5), une raideur apparait devant le terrne LU £ = 
V y • (a(x, y)V y U e ) qui permet de s’inspirer des methodes de developpement asymptotique largement 
utilisee en theorie cinetique pour construire des schemas numericiues multi-echelles (voir [13]). Nous 
decomposons en effet la solution du probleme augmente U £ sous la forme U £ (t, x, y) = F £ (t, x)+G £ (t, x, y), 
ou F £ est la projection ortlrogonale de U £ sur le noyau de L. Cette decomposition permet de reformuler 
de fagon equivalente le probleme augmente en un systeme d’equations micro-macro satisfait par F £ et 
G £ . Notons que ce systeme ne contient aucune approximation et reste exact pour toute valeur de e. 

Nous nous basons ensuite sur cette decomposition pour construire un schema numerique multi-eclrelles. 
Grace a une discretisation en temps semi-implicite (inspiree de [13]), un schema Asymptotic Preserving 
est alors obtenu. Ce schema necessite l’inversion au niveau numerique de l’operateur L (a x fixe), connne 
dans le cas de la resolution numerique du probleme homogeneise (voir [7]). 

Bien que la methodologie soit presentee ici en dimension quelconque, nous effectuons des tests numeriques 
uniquement en dimension 1. Le bon comportement du schema pour differents £ G (0,1] est mis en evidence 
en le comparant a une methode directe et a la solution du probleme homogeneise. L’analyse du cas rnulti- 
dimensionnel avec une methode d’elements finis est en cours d’etude. Cette note est une version abregee 
du travail plus detaille [6]. 


1. Introduction 

For T > 0 and a smooth bounded domain f l C K d , we consider the following class of parabolic problems 

d t u £ = V x • [a(x,x/s)W x u £ ] + t€(0,T), x G ft, (1) 

where u e (t = 0,x) = g(x) is a given initial condition in L 2 (fl), / £ L 2 (0,T, L 2 (fl)) is a given source term, 
and we take for simplicity homogeneous Dirichlet boundary conditions u e = 0 in (0, T) x di1. The tensor 
a(x,y ) £ K dxd is assumed symmetric, uniformly elliptic and bounded, and periodic with respect to the 
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variable y = x/e £ Y = (0, l) d . The homogenization analysis as e —> 0 of such a multiscale problem is 
well-known and can be done using a two-scale convergence analysis, see [4,5,3,10]. Problem (1) admits 
a unique solution u £ is the space L 2 (0, T; H/ (Cl)) which converges towards an asymptotic solution u° as 
s —> 0, 

u e —> u° strongly in L 2 (0, T; L 2 (Cl)), u £ u° weakly in L 2 (0, T; Hq(CI)), 
where u° £ L 2 (0, T; H/ (11)) solves an effective non-oscillatory problem of the same form as (1), 

-V y ■ [a(x,y)[V x u° + V„ui[] =0, x £ Cl,y € Y, 
d t u°-V x - j a(x, y)[V x u°+V v ui]dy =/, x £ Cl, 

with the same initial and boundary conditions for u° in (3) as for u e and involving the elliptic “cell 
problem” (2) with solution U\(t,x,-) £ H/ er (Y), periodic with zero average with respect to the second 
variable y. Taking advantage of the separation of micro and macro scales, numerical homogenization 
methods exploit the above homogenization result to compute efficiently the solution of (1) in the asymp¬ 
totic regime e —> 0. For instance, such an efficient method is the heterogeneous Multiscale Method (HMM) 
[7] (see also the review [1]) which relies on a coupling of micro and macro finite element methods ap¬ 
plied to (2)-(3) (see [2] in the context of parabolic multiscale problems (1)). Having a computational cost 
independent of the smallness of e, HMM permits to approximate not only the asymptotic solution u°, 
but also the oscillatory solution u e and its gradient S7 x u e for a fixed small e using the approximation 
u e (t,x) ~ u°(t,x) + eui(t, x, y — x/e) based on corrector techniques. This latter approximation of u e is 
consistent only for small values of e but not for £ close to 1. We also mention the “multiscale finite element 
method” (rnsFEM) [9] (see the review [8]), which permits to compute the oscillatory solution u £ by using 
an enriched finite element space, but where the computational cost grows as e —► 0. The aim of this 
paper is to introduce a micro-macro decomposition which permits to approximate u £ accurately for both 
regimes £ —> 0 or e ~ 1 at a cost independent of £, in the spirit of Asymptotic Preserving schemes [11]. 

Remark 1 The asymptotic parabolic problem (3) is non-stiff and can be written in the form dtu° = Du°+f 
where D<f> = V x ■ (a°V x <f>) and a 0 is the so-called homogenized tensor. In other words, the asymptotic 
problem has the same form as (1) with a £ replaced by a 0 . In dimension d = 1, the homogenized tensor is 
given by the harmonic average a°(x ) = (f y a(x,y)~ 1 dy)~ 1 . However, in multiple dimensions, there is no 
such a simple formula. The calculation of the asymptotic diffusion coefficient a°(x) at a point x of space 
then requires the resolution of an elliptic type problem like (2). We refer to the review [1]. 

2 . Exact micro-macro decomposition 

In this section, we introduce a micro-macro decomposition of the oscillatory solution u e of (1) which 
remains exact for all £ £ (0,1) and we study its behavior for e —> 0. We emphasize that our analysis is 
only formal, and a rigorous study of the approach is currently under investigation [6]. In the spirit of 
two-scale convergence analysis, we introduce a function U £ : (0, T) x Cl x Y —► R, periodic with respect 
to the third variable y £Y = (0, \) d , and such that U £ coincides with u e , the solution to (1), on the 
diagonal, i.e. 

U £ (t, x, x/e) = u £ (t, x), (4) 

and we obtain that U £ solves the following augmented problem 

d t U £ = \LU e + -BU £ +DU e + /. (5) 

£ z £ 

Here, we use the following notations for all functions </>, 


( 2 ) 

( 3 ) 
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L(f> = Vy • [a{x, y)S7 v (f)], D<j> = V x • [a(x, y)V x <f>\, B<f> = V x • [a(x, y)V y <t>] + V y • [a(x, y)V x (f>\. (6) 

We then choose appropriate initial and boundary conditions for (5), 

U £ (0 , x , y) = g(x) for x £ Q,y £ Y, U £ (t , x, y) = £ ^«i(t, ar, y) — ui(£, x, —)^ for t G (0, T), a; € dft, y G Y, 

(7) 

where u\ is given by (2). Consider now the linear operator L in (6) defined for a fixed x on H[ er (Y ), the 
space of periodic functions in H 1 (Y); this operator is self-adjoint with respect to the L 2 {Y) scalar product. 
Its kernel is the set of constant functions (with respect to y) and the L 2 orthogonal projector on this 
kernel is the average projection operator He] := J Y (f>{y)dy. Moreover, L is an isomorphism between the 
Hilbert space W per (Y ) = {([ G H] er (Y) ; n tj> = 0} and its dual (W per (Y)Y. Following [13] in the context 
of kinetic theory, we now perform an exact micro-macro decomposition of U e by setting F e = n U £ , 
G £ = (I — n )U £ , where I is the identity operator 

U £ (t,x,y) = F £ (t,x) + G £ (t,x,y). (8) 

Inserting this decomposition into (5), and applying respectively n and (/ — n) leads to the following 
model for the micro-macro decomposition of u e (t,x) = F e (t,x) + G £ (t,x,x/e), 

d t F e = x n bg £ + n df £ + n dg £ + /, (9) 

£ 

F £ (0, x) = g(x), x G F £ (t, x) = — eui(t, x, —), t E (0, T),x G 

£ 

d t G £ = \lg £ + 1 (/ - n )b[f £ + g £ ] + (i - n )d{f £ + g £ ), ( 10 ) 

£ z £ 

G £ (0,x,y) = 0,x G Q,y G Y, G £ (t, x, y) = eui(£, x, y), t G (0, T), x G dtt, y G Y, 

where u\ is given in (2), and we note that HBcj) = 0 for all <f> independent of y. 

We now study formally the asymptotic limits of G £ and F £ —>• F as e —> 0. Multiplying both sides of 
(10) by £, using d t G £ = 0(1) and setting formally £ —¥ 0, we deduce 

e-'G 6 -> -L-'BF. (11) 

Injecting (11) into (9) we deduce formally the asymptotic limit as £ —> 0, 

d t F = DF + /, with Dcj) := -BBL~ 1 B(j) + UDcj). (12) 

Comparing with the homogenized problem (2)-(3), we emphasize that the asymptotic limits as £ —> 0 
in (11) and (12) coincide respectively with u° and u\ in (2),(3). Precisely, we have F e —► F = u° and 
s~ 1 G £ —L~ 1 Bu° = Ui for £ —»• 0, and we see that the parabolic effective problem (3) is equivalent to 

(12). Therefore, the decomposition (4)-(8) of u £ (t,x ) into the non-oscillatory part F £ (t,x) and the oscil¬ 
latory part G £ (t,x,x/e) can indeed be interpreted as a generalization for £ G (0,1) of the approximation 
u £ ~ u° + EUi valid only for small values of e. 

Remark 2 Notice that using the simpler homogeneous boundary condition U(t,x,y) = 0 for t G (0,T),a;G 
dtt,y G Y (instead of (7 )) would yield an undesired boundary layer. Indeed, it is a classical issue for 
corrector techniques (see [10]) that the limit U\(x,y) appearing in (11) with y = x/e does not satisfy 
homogeneous Dirichlet boundary conditions for x G dft. 

3. Asymptotic Preserving numerical method 

The main goal of this section is to propose an Asymptotic Preserving numerical method for (9)-(10), 
used to approximate the solution u £ of (1), and with a computational cost independent of £ G (0,1). We 
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focus on the time discretization, considering a mesh of the time interval [0,T]: t n = nAf, with nSN and 
At the time step. Then, we denote by F n (resp. G n ) an approximation of F e (t n ) (resp. G £ (t n )). 

The stiffest term in (10) has to be considered implicit to ensure stability as £ —> 0 whereas all the other 
terms can be considered explicit. Then, a natural first order time discretization of (10) is (see [13]) 


G n+1 = I I~ 


-1 r 


G n + — (I -U)(B + eD){F n + G n ) 


Now, to recover the correct asymptotic behavior, the time discretization of (9) is 

F n+ 1 =F n + —HBG n+1 + AtHD(F n + G n+1 ) + /. 

£ 


(13) 


(14) 


To make explicitly appear the asymptotic model, we now propose a suitable transformation of the macro 
part (14) (following [12]). To do that, we consider the following Duhamel formula for (10) 


d t {e- tL/e2 G e ) = e~ tL/s2 


-(I- n ){BF e + BG £ ) + (I- n ){DF e + DG e ) 
£ 


Integrating between t n and t n+1 and performing some first order (in time) approximation leads to 


( J . n -\-1 


o A tL/e 2 qti 


i 

£ 



s)L/e2 ds(I- H)BF n 


1 




e (t n+1 s) L /e 2 ds(I - n )BG n + 



- s)L/e2 ds{I - n ){DF n + DG n ). 


(15) 


Now, our goal is to derive an approximation of G n+1 which will be inserted in the right hand side of (14). 
This means that, in the approximation of G n+1 , any term of the order of At can be neglected if needed. 
In this spirit, the time integrals in (15) are approximated as follows. The first one is calculated exactly to 

get J* e 0" +1 -s)£/e“_ — g 2 (l — e AtL / E2 )_L -1 whereas the second and the third ones are approximated 
using a midpoint formula f*„ + ~ s ') L / e ds ~ At, e AtL /( 2e ). Additionally, for non small e, up to terms 

of order At, e AtL ^ s can be replaced by e~ At ^ s ; for small e, both e AtL ! e and e~ At ! e go to zero. Therefore, 
g A tL/e ma y replaced by e~ At / £ . Let us emphasize that all the approximations are consistent in time 
with the continuous model. Finally, we get 


G 


n+1 


= e- At/£2 G n -e{ 1-e 


- At / £2 )L~ 1 (I-U)BF n +—e- At/ ^ £2) (I-U.)BG n +Ate- At/ ^ £2) (I-IV)(DF n +DG n ). 

£ 


This expression of G n+1 is injected in (14) and using —II^BL 1 (/ — U)BF n + TLDF n = DF n , we have 

F n+ 1 =F n + At(l - e - At/£2 )DF n + Ate~ At/£2 UDF n + — - —llBG n + AtHDG n+1 + At/. (16) 

£ 

The numerical scheme (13)-(16) enjoys the Asymptotic Preserving property: (i) for a fixed £ > 0, it is a 
first order approximation of (9)-(10); (ii) for a fixed At, (13)-(16) is uniformly stable with respect to £ 
and degenerates into a consistent discretization of the asymptotic model (12). 

Note that up to first order terms in At, an implicit time discretization for the term DF could also be 
considered (or alternatively an explicit stabilized method as in [2]), which enables to avoid the parabolic 
CFL condition in the asymptotic regime. Let us also remark that the numerical scheme (13)-(16) only 
requires the inversion of elliptic type operator in dimension d, which is also needed for the numerical 
resolution of the asymptotic model (12). Hence, the additional computation induced by our approach 
only lies in the numerical approximation of n, B , D. 
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4. Numerical results 


We consider a simpler case to illustrate the efficiency of our approach, considering (1) in dimension 
d = 1 together with the following diffusion coefficient A{x/e) = 1.1 + sin(a;/e) and a null right hand side 
f(x) = 0 with the initial condition g(x) = sin(27r;r), x £ [0,1]. We then compare three different approaches 
for approximating the solution u £ to (1): (i) a direct approach (’REF’) based on a first order explicit time 
integrator of (1) in which the mesh parameters are adapted to the smallness of £ (this will serve as a 
reference for comparison, obviously, another time integrator can be used, such as an implicit one); (ii) the 
exact micro-macro decomposition approach (’EMM’) based on (13)-(16) and (in) the asymptotic model 
(’HMM’) u e ~ u° + £U\ with u°,u\ given by (2)-(3). The spatial discretization (in x and y directions) is 
performed using a standard second order finite volume method. We denote by N x (resp. N y ) the number 
of points in the x (resp. y) direction, and Aa: = 1 / N x . Ay = 1 /N y are the mesh size in x and y directions. 
For REF, we choose At = 0.05A;r 2 where Aa; will be chosen small enough to capture the oscillations 
of size £. For EMM and HMM, we choose At = 0.2Aa; 2 . Once U £ (t,x,y) is computed on the discrete 
meshes, we need to interpolate it at y = x/e, y being a periodic variable; this is done using trigonometric 
interpolation which ensures spectral accuracy. As a diagnostic for EMM, we consider a reconstructed 
solution on a refined mesh by using the following linear interpolation for x £ [a;,;, x,+i], i = 0,..., N x — 1, 

u s (t n ,x) » {Xi+ ^ x) [F n (x i )+G n (x i , j)) + (F n (x m ) + G n (x i+1 , *)) , (17) 

where ( F n ,G n ) given by (16)-(13), which enables to recover the small-scale information. For HMM, we 
use the same reconstruction (17) with u° and eu\ at t = t n (given in (2)-(3)) instead of F n and G n , 
respectively. Hence, this enables to have an approximation of the EMM and HMM solutions as well as 
their spatial derivative on a refined mesh on which the reference solution has been obtained. 

In Figure 1, we consider the cases e = 1,0.1,0.01, respectively (from left to right). We plot at the final 
time T = 1 the solutions given by REF, EMM and HMM as functions of x £ (0,1) (first line), the error in 
u e (second line, plotting u RB f — «emm and Wref — Uhmm), and the error in the derivative d x u £ (third line, 
plotting clcitREF — c^Wemm and c^Wref — 9 x u H mm)- For the reference solution REF, we use N x = 1024 
for e = 1,0.1, and N x = 4096 for e = 0.01. Concerning EMM and HMM, we use in all computations the 
mesh parameters N x = 64 and N y = 16 for all e. We can observe that for arbitrary values of £, EMM is 
in a very good agreement with the refined REF solution and its derivative, for a fixed set of numerical 
parameters. However, the HMM solution is accurate only in the case £ = 0.01. In addition, as £ goes 
to zero, the computational cost for EMM is constant whereas the one of a direct method such as REF 
increases as the meshes need to be refined. 
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(from left to right). Second line (errors in u e ): uref - ^emm (red line) and uref _ ^hmm (blue line) for e = 1, 0.1, 0.01 (from 

left to right). Third line (errors in d x u £ ): ^uref (red line) and ^uref — 9 s uhmm (blue line) for e = 1, 0.1, 0.01 

(from left to right). The EMM and HMM solutions and their derivatives are computed using the interpolation (17). 
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